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We present a first-principles study of ferroelectric domain walls (FE-DWs) in multiferroic BiFe03 
(BFO), a material in which the FE order parameter coexists with anti-ferrodistortive (AFD) modes 
involving rotations of the 06 octahedra. We find that the energetics of the DWs are dominated 
by the capability of the domains to match their C>6 octahedra rotation patterns at the plane of 
the wall, so that the distortion of the oxygen groups is minimized. Our results thus indicate that, 
in essence, it is the discontinuity in the AFD order parameter, and not the change in the electric 
polarization, what decides which crystallographic planes are most likely to host BFO's FE-DWs. 
Such a result clearly suggests that the 06 rotational patterns play a primary role in the FE phase 
of this compound, in contrast with the usual (implicit) assumption that they are subordinated to 
the FE order parameter. Our calculations show that, for the most favorable cases in BFO, the 
DW energy amounts to a several tens of mj/m 2 , which is higher than what was computed for 
other ferroelectric perovskites with no 06 rotations. Interestingly, we find that the structure of 
BFO at the most stable DWs resembles the atomic arrangements that are characteristic of low-lying 
(meta)stable phases of the material. Further, we argue that our results for the DWs of bulk BFO 
are related with the nanoscale-twinned structures that Prosandeev et al. [Adv. Funct. Mats. (2012), 
doi: 10. 1002/adfm. 201201467] have recently predicted to occur in this compound, and suggest that 
BFO can be viewed as a polytypic material. Our work thus contributes to shape a coherent picture 
of the structural variants that BFO can present and the way in which they are related. 

PACS numbers: 77.84.-s, 77.80.Dj, 75.85.+t, 71. 15. Mb, 61.50.Ah 



I. INTRODUCTION 

Ferroelectric (FE) crystals are insulators displaying a 
macroscopic polarization that can be switched by the ac- 
tion of an electric field. The ground state of a bulk fer- 
roelectric, where a unit cell is periodically repeated in 
space, is a perfectly ordered structure where the atomic 
patterns that give rise to the polarization are the same 
in every part of the crystal. However, real ferroelectric 
crystals need to minimize the electrostatic energy aris- 
ing from their finite size (i.e., the electric field outside a 
ferroelectric sample must rapidly decay to zero), and as 
a result the formation of regions where the polarization 
points in different directions will occur. Like in the case 
of other ferroic materials such as ferromagnets or ferroe- 
lastics, these regions are called domains, and the area in 
between domains is called domain wall (DW). 

The presence of DWs in a material can cause important 
changes to its macroscopic behavior P T hus , for example, 
recent findings about the conducting^! and photovoltaic 4 
properties of ferroelectric DWs (FE-DWs) have fostered 
the interest in using them as functional parts in devices 
in nanoelectronics!^ The vast majority of the recent dis- 
coveries on FE-DWs have featured multiferroic BiFe03 
(bismuth ferrite or BFO), a material that has concen- 
trated many experimental efforts during the past decade 
because of the promise it holds for room-temperature ap- 
plications based on magnetoelectric effectsPFor all these 



reasons, currently there is a strong interest in improving 
our atomistic understanding of FE-DWs, in particular in 
the case of BFO. 

While the width of magnetic DWs is in the range of 
10 to 100 nm, FE-DWs are believed to be much thinner. 
The atomic structure of ferroelectric walls can be studied 
experimentally using techniques such as high-resolution 
transmission electron microscopy, which permits a spatial 
definition of up to around 0.5 A. However, the atomic dis- 
placements that characterize FE-DWs are of the order of 
0.02 A, and therefore direct imaging and interpretation 
of the structure of these walls is still a challenge.^ First- 
principles calculations have also been used to obtain in- 
formation about FE-DWs. Some studies^!3 were carried 
out on ferroelectric perovskites where the polarization 
is the only structural order parameter, such as BaTiOa 
and PbTi0 3 . They report that FE-DWs in these ma- 
terials are atomically thin. Other recent first-principles 
studies2H2i have focused on multiferroic BFO. The width 
computed for BFO's DWs is at most 1 nm, and some of 
them are reported to show a chiral pattern that is abstent 
in former studies. 

In the rhombohedral phase of BFO that is stable at 
ambient conditions (i?3c space group), the Oq octahedra 
of oxygens that surround the Fe ions are rotated in an- 
tiphase about the polarization axis. These rotations con- 
stitute a second structural order parameter (OP) in addi- 
tion to the polarization; such an 06-rotation related OP 
is usually referred to as anti-ferrodistortive (AFD). AFD 
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modes constitute the most usual structural instability in 
the family of perovskite oxides, and have been studied in 
a systematic and extensive way starting with the classic 
works of Glazei^ and others. Interestingly, AFD and 
FE distortions often compete in perovskite oxides, and it 
is rare to see them occur simultaneously in a material.^ 
In fact, as far as we know, BFO is the only perovskite 
oxide that presents both OPs in its equilibrium phase 
at ambient conditions, which makes it a unique case for 
study. 

Whenever two OPs occur simultaneously in a given 
phase of a material, we tend to assign the lables primary 
and secondary to them. Such a classification is useful 
to rationalize the structural phase transition and domain 
variants occurring in the crystal: The primary OP is the 
one that drives the phase transition and it would occur 
even in absence of the secondary OP; the primary OP 
also defines the number and properties of the structural 
domains that can form. In contrast, the occurrence of 
the secondary OP relies on the presence of the primary 
distortion; the secondary OP may appear as a result of 
the symmetry breaking caused by the primary OP (its 
character would thus be improper) or because of a desta- 
bilizing coupling with the primary OP (the secondary 
order would be triggered in that case). Interestingly, as 
we will see, such a distinction between primary and sec- 
ondary OPs cannot be made in BFO's case. 

Most of the literature on BFO focuses on the investi- 
gation of the FE polarization, which is effectively treated 
as the single structural OP of the compound. This is a 
rather natural approach: The spontaneous polarization, 
with its associated electrostatic energy, is the leading 
driving force for the formation of domains; additionally, 
we can easily act on it by applying an external field. As 
regards the AFD distortions, the oxygen octahedra are 
known to rotate specifically about the polarization axis, 
and are much more difficult to manipulate and charac- 
terize experimentally; hence, they are usually treated as 
if they played a secondary role. However, these AFD 
modes do not occur as a result of a polarization-induced 
symmetry breaking; further, first-principles calculations 
show that they constitute a genuine structural instabil- 
ity of the ideal cubic perovskite phase of BFO, indepen- 
dently of the occurrence of the polarization. Hence, it 
is clear that the AFD distortions do not fit the defi- 
nition of a secondary OP given above. Moreover, the 
rhombohedral phase of BFO does not occur as a conse- 
quence of a displacive phase transition characterized by 
the condensation of the FE soft mode; rather, BFO's i?3c 
phase undergoes a strongly discontinuous transformation 
upon heating, and the resulting high-temperature struc- 
ture is probably characterized by oxygen-octahedra rota- 
tions more than it is by polar distortions. Hence, from 
this perspective either, there are no clear arguments to 
distinguish between primary and secondary OPs in BFO. 

Detailed studies based on density functional theory 
(e.g., see Table I of Ref. [THJ) allow us to make this 
discussion more quantitative. It has been shown that 



the FE distortion of the ideal perovskite phase (which 
would lower the symmetry from the cubic Pm3m space 
group to rhombohedral R3m) reduces BFO's energy by 
about 750 meV/f.u. (the precise result depends on the 
density functional employed in the calculations^). In 
contrast, a pure AFD distortion (leading to the i?3c 
space group) reduces the energy of the cubic phase by 
about 680 meV/f.u. The combined FE and AFD dis- 
tortions render the Ric phase of BFO, which lies about 
925 meV/f.u. below the cubic structure. We can thus 
conclude that the FE and AFD instabilities compete: If 
they did not interact at all, the energy of the phase com- 
bining both should be about 1430 meV/f.u. below the 
cubic structure (with 1430 = 750+680); the significantly 
smaller energy difference obtained from the simulations 
(925 meV/f.u.) clearly denotes a competition. 

Hence, the simultaneous occurrence of FE and AFD 
distortions in BFO's R3c phase has to be attributed to 
the fact that, individually considered, both constitute 
very strong instabilities of the cubic phase of the ma- 
terial. Such instabilities clearly compete, not cooperate; 
yet, their simultaneous occurrence renders a net reduc- 
tion of BFO's energy, thus overcoming the adverse effect 
of the competition. Let us note that, because the po- 
larization and rotation axes coincide, one might be led 
to believe that the two OPs cooperate in some way, but 
such an interpretation would be wrong; rather, the coin- 
cidence of the two axes defines the conditions for which 
the FE-AFD competition is weakest. 

In view of these facts, one should probably say that 
BFO presents two primary order parameters that need 
to be considered on equal footing in a theoretical discus- 
sion of this material. Bearing this in mind, we have con- 
ducted a first-principles study to reexamine the FE-DWs 
occurring in the R3c phase of bulk BFO. Interestingly, 
we find that the DW energetics is essentially determined 
by the type of discontinuity affecting the AFD patterns, 
and not so much by the change in the FE polarization 
accross the wall. Our observations have a rather natu- 
ral explanation in the context of recent discoveries high- 
lighting what we may call BFO's polymorphism, as we are 
able to relate the lowest-energy DWs found with the low- 
lying meta-stable phases that this compound is known to 
present. More precisely, the atomic arrangement at the 
most stable FE-DWs resembles the structure of a dif- 
ferent BFO polymorph. Thus, our results suggest that 
multi-domain configurations in BFO display a peculiar 
form of polytypism. The connection of our work with the 
novel nanoscale-twinned structures^ predicted to occur 
in BFO under various contidions (high temperature, high 
pressure, appropriate chemical substitution) are also dis- 
cussed, offering a coherent picture of structural diversity 
and energetics in this material. 

The paper is organized as follows. In Section II 
we discuss the details of the configurations we studied, 
introducing a rigourous classification of the FE/AFD 
DWs that can occur in BFO; we also describe the first- 
principles methodology used. In Section III we present 
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(b) 

FIG. 1. (Color online.) (a) Pseudocubic cell of BFO showing 
its three main axes (thin arrows) and the polarization along 
[lll] P c (thick arrow), (b) Detail of the cell, where two octa- 
hedra have been removed to allow visualization of Fe cations 
and Fe-O bonds. 

and discuss our results, revealing the structural details 
that determine the DW energetics. Other aspects of 
the DWs, as e.g. their electronic structure, are also ad- 
dressed, and the connection of our results with existing 
experimental information is discussed. In Section IV we 
summarize the work and present our conclusions. 

II. METHODS 

A. Definition of the combined FE/AFD-DWs 

BFO crystallizes in rhombohedral space group i?3c, 
with a primitive unit cell that contains 10 atoms. Its 
rhombohedral lattice parameters are a r h = 5.6343 A and 
ct r h = 59.348° It is easier to visualize its structure us- 
ing the pseudocubic 40-atom unit cell shown in Fig. [T] 
The Bi and Fe cations are displaced from their high sym- 
metry positions along the [lll] pc pseudocubic axis, and 
the material develops a polarization in this direction. At 
the same time, the 06 octahedra surrounding Fe cations 
rotate about this [lll] pc axis in antiphase, which cor- 
responds to the so-called a~a~a~ pattern in Glazer's 
notation.^ 

Now, the cubic perovskite structure presents 8 equiv- 
alent directions of the (lll) pc type. In a given domain, 
the electric polarization P can point along any of these 
8 directions; further, once a direction for P is chosen, 



there are two equivalent ways (related by a phase shift) 
in which the Oq octahedra can rotate. Hence, in total we 
can have 16 different FE+AFD domains in the rhombo- 
hedral phase of BFO. 

We now introduce a notation that will help us to de- 
scribe the relationships between these domains. Let u>i 
be the three-dimensional pseudovector quantifying the 
rotations of the Oq octahedron at cell z; the components 
U!i a , with a — x,y, and z, pertain to individual rota- 
tions about each of the pseudocubic principal axes of 
the perovskite structure, which are sketched in Fig. [T] 
(From now on, all vectors and planes will be referred 
to the pseudocubic axes unless we mention otherwise!^) 
As mentioned earlier, BFO's Og rotations around (111) 
arc modulated in antiphase when we move from one 
cell to any of its first-nearest-neighboring cells. Such 
a modulation corresponds to the R g-point of the Bril- 
louin zone of the ideal cubic 5-atom perovskite cell; for 
the particular case of rotations around the [111] axis, 
the pattern can be mathematically expressed as Wj = 
uj (— 1)""+™^+™" (l ; l ; l)j where the n ia variables are in- 
tegers defining the location of cell i in the lattice. 

Let us consider how a FE-DW can affect this ideal pat- 
tern. To do that, let us choose a cell i that is well inside 
the first domain; let us assume that u;j = wo(l, 1, 1) and 
that the polarization in this first domain is also along 
[111], so that P 1 = P (l, 1, !)• Then, let us pick a cell 
i' that is well inside the second domain, and which can 
be reached from cell i by advancing an even number of 
cells along the pseudo-cubic direction a. Thus, we have 
nvp — n ifi = ^nd a /3, n being an integer. Obviously, the 
two domains separated by the DW are related by sym- 
metry, which implies that uij and a>j< must be symmetry- 
related as well. We can have the following possibilities: 

(1) Wj' = o?j, which coincides with what we would see 
in the monodomain situation, i.e., the pattern of 06 ro- 
tations is essentially unaltered by the FE-DW. We will 
denote this case as a "AFD-DW of type 0". Note that, 
since in BFO the 06-rotation axis coincides with the di- 
rection of the spontaneous polarization, we can conclude 
that the polarization of the second domain must be ei- 
ther P n = P (l,l,l) or P n = P (-l,-l,-l); the for- 
mer possibility implies there is no FE-DW ("FE-DW of 
type 0") and the latter corresponds to a 180° FE-DW 
( "FE-DW of type 3" , as the three polarization compo- 
nents change sign). We have thus identified the first two 
possible types of FE/AFD-DWs that can occur in BFO, 
and which we will denote, respectively, by "0/0" (i.e., FE- 
DW of type and AFD-DW of type 0, which means that 
there is no discontinuity at all) and "3/0" (i.e., FE-DW 
of type 3 and AFD-DW of type 0, which means that the 
discontinuity only affects the polarization, which rotates 
by 180° when moving from domain I to domain II). 

(2) One rotation component changes sign with respect 
to the ideal situation. For the sake of the discussion, 
let us assume that the affected direction is z, so that 
u>j/ = Wrj(l, 1) — 1). This case, which we will denote as 
"AFD-DW of type 1" , does imply a discontinuity of the 
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a~a~a~ rotational pattern. Indeed, we must have some 
sort of phase boundary for the Oq rotations about the 
z axis as we cross the DW. As we will see later, our 
simulations indicate that such a phase boundary is the 
most trivial one that we could imagine: The z-rotations 
of neighboring cells accross the DW occur in-phase, thus 
breaking the anti-phase modulation sequence of the ideal 
pattern. Here again, we can deduce that the polarization 
of the second domain must be either P n = Po(l, 1, — 1) 
or P n = P (-l, -1, 1) implying that the FE-DW must 
be either of "type 1" (the 71° FE-DW in which only 
one polarization component changes sign) or of "type 2" 
(the 109° FE-DW with two components changing sign), 
respectively. We have thus identified two new types of 
FE/AFD-DWs, which we denote by "1/1" and' "2/1", 
respectively. 

(3) Two rotation components change sign, so that e.g. 
uv = wo(l, -1, -!)• This will be our "AFD-DW of type 
2", and the associated FE-DW can be either of "type 
2", with P 11 = P Q (1,-1,-1), or of "type 1", with P n = 
Po(-l, 1, 1). We thus have the combined FE/AFD-DWs 
of types "2/2" and "1/2". 

(4) u>j/ = — Ui, corresponding to a "AFD-DW of type 
3". In this case, the FE-DW can either be of "type 
3", with P n = P (-l,-l,-l), or of "type 0" , with 
P n = P (l, 1, 1). The combined FE/AFD-DWs are thus 
denoted by "3/3" and "0/3", respectively. Obviously, in 
the latter case there is no FE-DW, and the discontinuity 
only affects the AFD pattern. 



B. Domain walls investigated 

In this article we will consider DWs occurring at planes 
of low Miller indices, namely, (100) and (110), which are 
the ones that were discussed by previous authors Our 
initial task is to reduce the number of possible combi- 
nations of domains that we have to study. First, for a 
given DW many of these configurations will be equiva- 
lent by symmetry. Second, we limit ourselves to neutral 
DWs, i.e., those in which there is no discontinuity in 
the polarization component across the wall. Note that 
such a discontinuity would give raise to bound charges 
at the DW, whi ch wo uld lead to a large unfavorable elec- 
trostatic energy.E2122l Hence, the electrostatically allowed 
cases are those in which (P 1 — P n ) • (1,0,0) = and 
(pi _ pii) . (1,1,0) = 0, respectively, for the (100) and 
(110) planes. 

Lastly, we could also limit ourselves to dom ains t hat 
fulfill the condition of mechanical compatibility)^^ To 
understand the origin of this condition, note that the 
ground state phase of ferroelectric perovskites involves a 
distortion of the cubic cell, i.e., an homogeneous strain. 
In the presence of a DW, the distortion of unit cells at 
opposite sides of the DW may or may not coincide. If 
such strains coincide, we have mechanical compatibility; 
if not, an elastic energy penalty will occur. In the case of 
BFO, a possible mechanical mismatch at the DW comes 



TABLE I. Configurations studied. Columns: (1) polarization 
direction in domain I, (2) DW plane, (3) polarization direction 
in domain II, (4) FE/AFD DW-type according to polarization 
change and octahedra rotation pattern change (see text), and 
(5) whether the mechanical compatibility condition is met at 
the DW. 



P 1 


DW 


pii 


Type 


Mech.? 


[111] 


(100) 


[111] 


0/0 


YES 


[111] 


(100) 


[111] 


0/3 


YES 


[111] 


(100) 


[111] 


1/1 


NO 


[111] 


(100) 


[111] 


1/2 


NO 


[111] 


(100) 


[111] 


2/1 


YES 


[111] 


(100) 


[111] 


2/2 


YES 


[111] 


(110) 


[111] 


0/0 


YES 


[111] 


(110) 


[111] 


0/3 


YES 


[111] 


(110) 


[111] 


1/1 


YES 


[111] 


(110) 


[111] 


1/2 


YES 


[111] 


(110) 


[111] 


0/0 


YES 


[111] 


(110) 


[111] 


0/3 


YES 


[111] 


(110) 


[111] 


1/1 


NO 


[111] 


(110) 


[111] 


1/2 


NO 


[111] 


(110) 


[111] 


2/1 


NO 


[111] 


(110) 


[111] 


2/2 


NO 


[111] 


(110) 


[111] 


3/0 


YES 


[111] 


(110) 


[111] 


3/3 


YES 



from the relatively small shear strains that change the 
pseudocubic angles from 90° to about 89.2°. As a result, 
in BFO there are pairs of domains that, although do not 
exactly fulfill the mechanical compatibility condition, are 
very close to doing so. We thus decided to consider them 
for study. 

Once the criteria described in the previous paragraphs 
are applied, it turns out that there are 6 possible inequiv- 
alent configurations for (100) DWs, all of them with P 1 = 
Po(l, 1, 1). Concerning (110) domains, there are two in- 
equivalent ways to fix domain I: using P 1 = Pq(1, 1, 1), or 
using P 1 = Pq(1,— 1,1), which implies a domain where 
the polarization is parallel to the DW. For these two al- 
ternatives we find, respectively, four and eight inequiv- 
alent DW configurations. In total, we have 18 different 
DWs that we considered in this study. Table [T] lists them, 
indicating their FE/AFD-DW type and whether the me- 
chanical matching condition is fulfilled or not. 



C. First-principles methods 

Our calculations are based on density-functional the- 
ory (DFT). 2 " 3 Most of them were done using the local- 
density approximation (LDA)(23H2H f or the exchange 
and correlation functional, although for comparison we 
repeated some calculations using two flavors of the 
generalized-gradient approximation: the functional of 
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Perdew, Burke, and Ernzerho^ (PBE), and its adap- 
tation to solids 28 (PBEsol). To obtain a better descrip- 
tion of iron's 3d orbitals, we added a Hubbard-[/ term 
to the energy of the system, following the prescription of 
Dudarev and coworkers.^ 

We used the implementation of this formalism in the 
Siesta codeP^ Siesta uses a basis set of localized nu- 
merical orbitals, which makes it computationally very ef- 
ficient. In our case, we included in our basis set s (DZ), p 
(SZ), and d (DZ) orbitals centered at Fe; s (DZ), p (DZ), 
d (SZ), and / (SZ) orbitals centered at Bi; and s (DZ), 
p (DZ), and d (SZ) orbitals centered at O. The SZ and 
DZ terminology means that we used either one (singlc- 
C) or two (double-^) radial functions for each of the 
angular forms, respectively.^ We used norm- conserving 
Troullier-Martins pseudopotentials^ for a efficient treat- 
ment of the interaction between ion cores and valence 
electrons; the pseudopotentials were generated in the fol- 
lowing configurations: Fe's 4s 2 3d 6 , Bi's 6s 2 6p 3 , and O's 
2s 2 2p 4 ; partial core corrections were included in all cases. 
We used a real-space grid to represent the electronic den- 
sity, the Hartree and exchange-correlation potentials, and 
the matrix elements between basis orbitals; this grid is 
equivalent to one generated by a plane- wave basis cut off 
at 300 Ry. The integrations in reciprocal space were done 
using a Monkhorst-Pack grid equivalent to a 4 x 4 x 4 one 
for a 5-atom perovskite unit cell. 

We also used the Vasp code 3 - 2 ^ for our calculations. 
Vasp uses plane-waves as basis sets, and therefore al- 
lows for a more systematic numerical treatment of the 
equations to be solved at the cost of more computational 
time. We used the projector-augmented wave methocP^ 
to represent the interaction of electrons with ionic cores 
solving for the following electrons: Fe's 3s, 3p, 3d, and 
4s; Bi's 5d, 6s, and 6p; and O's 2s and 2p. We generated 
a plane-wave basis set using a kinetic energy cutoff of 
500 eV. The integrations in reciprocal space were done 
using the same k-point grids as in Siesta. 

When using Vasp we set the Hubbard-JJ equal to 4 eV, 
a value determined by requesting the computed magnetic 
interactions to be in quantitative agreement with those 
obtained from calculations with hybrid functional 
U ~ 4 eV has become a frequent choice in first-principles 
studies of BFO with this same kind of formalism, as it 
leads to qualitatively and semi-quantitatively correct re- 
sults for all the properties investigated so far. Due mainly 
to the differences in the treatment of the ion cores, with 
Siesta we need to use U = 2 eV to reproduce the Vasp 
results for the ground state of BiFeC>3; thus, we used this 
value for the Siesta calculations. 

Unless mentioned otherwise, the results in our tables 
and figures are the ones obtained with Vasp. 



TABLE II. Energetics of the configurations studied with sim- 
ulation cells containing 80 atoms. Column 1: domain type 
described following Table [T] Column 2: DW energy according 
to our simulations; asterisks indicate configurations where the 
mechanical matching condition is not fulfilled. Column 3: la- 
bel according to Oq matching pattern (see text). 



DW Type E ow (mj/m 2 ) Matching 



[1111 

L J 


(100) 


[1111 


0/0 







[1111 

[X XXj 


(100) 


[1111 

L J 


0/3 


227 


c 


[1111 
[XXXJ 


(100) 


[1111 

[XXXJ 


1/1 


151* 


B 


[1111 

L J 


(100) 


[1111 

L J 


1/2 

X/ £j 


147* 


B 


[1111 

L J 


(100) 


[1111 

L J 


2/1 

Aj X 


62 


A 


[1111 

L J 


(100) 


[1111 

L J 


2/2 


319 


c 


[in] 


(110) 


[111] 


0/0 







[in] 


(110) 


[111] 


0/3 


254 


c 


[in] 


(110) 


[111] 


1/1 


152 


B 


[in] 


(110) 


[111] 


1/2 


178 


B 


[Hi] 


(110) 


[111] 


0/0 







[Hi] 


(110) 


[111] 


0/3 


293 


C 


[Hi] 


(110) 


[111] 


1/1 


142* 


B 


[Hi] 


(110) 


[111] 


1/2 


196* 


B 


[Hi] 


(110) 


[ill] 


2/1 


150* 


B 


[Hi] 


(110) 


[ill] 


2/2 


244* 


B 


[Hi] 


(110) 


[111] 


3/0 


74 


A 


[Hi] 


(110) 


[ill] 


3/3 


255 


C 



III. RESULTS AND DISCUSSION 

A. Energetics and DW structure 

As a first step in characterizing our DWs we computed 
the energy needed to create them. In order to do so we 
designed an 80-atom periodic simulation cell for each con- 
figuration listed in Table [TJ containing two domains of 40 
atoms and two DWs. We then wanted to find the min- 
imum energy configuration for each case without mak- 
ing any prior assumption on the location of the DW in 
terms of parallel ion planes, nor on the specific atomic 
arrangement occurring at the DW. Thus, for each DW 
type considered, we started with some reasonable atomic 
configuration corresponding to the two domains involved, 
and performed a short molecular dynamics simulation, 
starting from random velocities, using Siesta. Our tests 
indicated that this procedure leads to an unrestricted 
and efficient exploration of the energy surface, allowing 
the system to find the lowest-energy DW configuration 
in a reliable and reproducible way. We then relaxed the 
lattice vectors and atomic positions until the stresses and 
forces were small. The configurations thus obtained were 
further relaxed using Vasp until the forces were below 
0.01 eV/A and the stresses were below 0.1 kbar. At all 
times the magnetic moments of the Fe ions kept an anti- 
ferromagnetic arrangement of type G (closest neighbors 
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have antiparallel spins); this is the configuration that 
BFO shows experimentally at a local scale at ambient 
conditions. Table [TT] includes all the DW energies that 
resulted after this process, computed as 



E — Eq 
2S 



(1) 



where E is the energy of the DW configuration, Eq is the 
energy of bulk BFO (computed for the same supercell 
for better numerical accuracy), S is the surface area of 
the cell face parallel to the DW, and we divide by two 
because each cell contains two DWs. 

Before discussing our results, let us note that the £?dw 
values that we obtain are significantly lower that those re- 
ported in the previous studies of Seidel et aZP and Lubk 
et aZP^ A comparison is shown in Table III where we 



include our results for the lowest-energy DWs that in- 
volve polarization rotations by 71°, 109°, and 180°, re- 
spectively. For consistency with the previous studies, all 
calculations reported in this Table were done with cells 
that contained 120 atoms, and we have restricted our- 
selves to cases satisfying the mechanical compatibility 
conditions. (Our DW energies obtained using 120-atom 
cells are very similar to the ones we calculated with 80- 
atom cells, implying that the results of Table [TT] are well 
converged in this regard.) There are some minor dif- 
ferences between the Vasp calculations done by the au- 
thors of Ref s . [21 and [T2l and ours; yet, we have tested that 
even if we modify our numerical approximations to make 
them as similar as possible to theirs, our results change 
very little. Moreover, Table |III| also includes the results 
we obtaind with Siesta, which are very similar to the 
values found with VASP. In addition, let us note that 
similarly low DW energies have been obtained by Wei et 
a ]|35] m an investigation of FE-DWs in BFO thin films un- 
der various epitaxial strain conditions. What can explain 
the discrepancies with the previous published literature? 
As some of us have recently shownpSJ the energy surface 
of BFO presents a very large number of local minima, 
which seems to be mainly a consequence of the many 
bonding complexes allowed by Bi's chemical versatility. 
We reached exactly the same conclusions in the present 
investigation of BFO's DWs, as we found that many pos- 
sible DW configurations can render a local minimum of 
the energy surface. Naturally, the physical relevance of 
high-lying local energy minima is questionable, and we 
found ourselves forced to implement a careful search for 
the global energy minimum for each DW type. Indeed, 
by thermally agitating and then quenching our configura- 
tions during the optimization process, we have had access 
to minima with (much) lower energy than the ones previ- 
ously reported. Hence, our results stand a better chance 
of being physically relevant. 

We found that the [111] (100) [111] 2/1 DW configura- 
tion has the lowest DW energy. It occurs when two polar- 
ization components and one rotation component change 
sign at the DW (this is a 109° DW from the perspective 
of the electric polarization). The relaxed atomic struc- 



TABLE III. Energetics of the configurations with minimum 
energy for each type of FE-DW that fulfills the mechani- 
cal matching condition (simulation cells with 120 atoms): 
[111](100)[111] 2/1 (109°), [111](110)[111] 3/0 (180°), and 
[111](110)[111] 1/1 (71°). 







£dw (mJ/m 2 ) 






This work 


This work 


Lubk et sP 


DW Type 


(Siesta) 


(Vasp) 


(Vasp) 


109° 


63 


62 


205 


180° 


86 


82 


829 


71° 


197 


167 


363 



ture that we obtained is shown in Fig. [2^ a); here one 
can already appreciate that the DW discontinuity is very 
sharp, and the Og octahedra are not strongly distorted at 
or near the DW plane. Figures [3]ja) and|3jb) contain the 
values of the displacements of Bi and Fe cations along 
the pseudocubic axes, relative to the center of mass of 
the corresponding surrounding O12 dodecahedron and 0$ 
octahedron, respectively. Such distortions are the origin 
of the local polar dipoles giving raise to BFO's macro- 
scopic polarization, and thus allow us to monitor the P- 
change at the DW. The Fe/Bi displacements are plotted 
as a function of the distance between the corresponding 
cations and the center of the DW, which is assumed to 
be located in between the two planes of Fe ions whose en- 
vironment is most different from the one they see in bulk 
BFO. We can thus see that these polar displacements are 
almost identical to the ones in bulk, except for small dis- 
tortions in a very narrow area around the DW center. 
The 06 rotations are quantified in Figs, ffic) and[3^d); 
they are the three components of (— iy i i^+ n iy+ n "( J j ij m 
the notation of Section II. A, with the two graphs contain- 
ing information for the two inequivalent lines of neighbor- 
ing 06 octahedra that run towards the DW. The 06 rota- 
tion patterns are also very similar to the bulk ones, except 
for the discontinuity of one AFD component at the DW. 
Concerning Fe-0 distances, Fig. |3^e) shows again minor 
distortions confined to the DW. Finally, in Fig. [3jf ) we 
quantify the distortions of the 12 0-0 bonds that form 
the edges of each Oe octahedron; we do so by computing 
the root mean square of the differences between each of 
the 0-0 bond lengths in each 06 octahedron, dj, and 
their average length d, 



\ 



(2) 



This is zero for a regular octahedron, but takes a value of 
about 0.1 in bulk BFO, as the polar displacements cause 
variations in the lengths of the 06 edges. We see that 
this measure shows very small differences in 06 octahedra 
distortions with respect to the bulk ones. In all, the result 
of the detailed analysis of the atomic structure points to 
a narrow DW whose thickness is around 1 nm. 
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FIG. 2. (Color online.) Atomic structures obtained upon relaxation of 120-atom unit cells: (a) [lll](100)[lIT] 2/1 DW 
configuration; (b) [ill] (110) [111] 3/0 configuration; (c) [111] (110) [111] 1/1 configuration. The vertical lines indicate the 
approximate position of the center of the DWs. 



Our finding that the most stable DW configuration is 
precisely [111] (100) [111] 2/1 might seem somewhat ar- 
bitrary; on the contrary, this result leads to very sug- 
gestive conclusions. In this DW, the AFD component 
that changes sign is the one perpendicular to the wall; 
thus, if we focus on the two planes of octahedra next 
to the DW, we see that they display a a + b~b~ Glazer 
rotation pattern, i.e., the Oq tilts occur in phase (+ su- 
perscript) about the [100] direction perpendicular to the 
DW and in antiphase (— superscript) about the [010] 
and [001] directions within the DW plane. Further, the 
discontinuity in the FE polarization - which is mainly 
captured by the Bi displacements shown in Fig. [3]ja) - 
occurs in the plane of the DW, which results in an anti- 
polar pattern of Bi displacements around the DW cen- 
ter. Interestingly, such structural features are exactly the 
ones characterizing the Pnma phase of bulk BFO, which 
is experimentally known to occu r at h igh temperatures^ 
and under hydrostatic pressure) 37 * 38 ! and which is rela- 
tively close in energy to the Ric ground state according 
to previous first-principles calculations.^ In fact, BFO's 
Pnma phase would match the structure of our DW al- 
most identically if we added to it a polar distortion along 
the direction of the in-phase oxygen-octahedra rotations. 



Interestingly, first-principles theory suggests that such 
a ferroelectrically-distorted Pnma phase, which would 
present the polar space group Pna 2x, m ay be a low- 
lying meta-stable polymorph of BFO ESEH although such 
a structure has not been observed experimentally yet. 
Hence, our results suggest that the lowest-energy DW 
found in BFO owes its stability to the fact that it can 
adopt an atomic arrangement that mimicks a low-energy 
polymorph of the bulk material. 

Interestingly, our results for this DW are also reminis- 
cent of the so-called nanoscale-twinned phases that have 
been predicted to occur in BFO under various conditions 
(e.g., high temperature, hydrostatic pressure, chemical 
substitution of Bi by rare-earth cations).^ In fact, the 
novel phases described by Prosandeev et ai. in Ref. [THl 
which are claimed to act as a structural bridge between 
the i?3c and Pnma structures that appear in BFO's 
phase diagram, can be viewed as a sequence, along the 
[100] direction, of DWs of the type that we have just dis- 
cussed. It is thus conceivable that, as we heat up BFO's 
R3c phase, the Pnma structure will nucleate at DWs 
as the ones just described, giving raise to intermediate 
bridging phases of a polytypic nature. 

The second DW configuration that we found to display 
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FIG. 3. (Color online.) Values of magnitudes that character- 
ize the atomic displacements in the [111](100)[111] 2/1 DW 
configuration (polarizations at an angle of around 109°). The 
discontinuous lines show the values computed for bulk BFO. 
The continuous vertical lines mark the approximate position 
of the center of the DW. We describe in the text those mag- 
nitudes that are not obvious from the vertical axis labels. 



a similarly low DW energy is the [111] (110) [111] 3/0 DW, 
which is sketched in Fig. ^b). This is a 180° FE-DW 
where the cation displacements go in opposite directions 
at opposite sides of the wall, as shown in Figs. |4^a) and 
4|b). As regards the network of Og octahedra, there is no 
AFD discontinuity involved in this boundary; yet, as we 
can see in Figs. Qc) and|4^d), because the lattice needs 
to accommodate the large change in the polar cation dis- 
placements, the 06 rotation angles at the DW deviate 
noticeably from the bulk-like values. From Fig. |4^e) we 
notice that there are some Fe-0 bond lenght variations 
accompanying these distortions. However, the 0-0 bond 
length dispersion is very similar to the bulk one, as il- 
lustrated in Fig. |4|f). Let us also note that Lubk et 
al. computed a very high energy of 829 mJ/m 2 for their 
most stable 180° FE-DW which, as in our case, involves 
no AFD phase boundary; they also reported that the po- 
larization rotates in two steps from one domain to the 
other, forming a sort of chiral pattern. Our finding of 
a sharp 180° FE-DW with an energy of 82 mJ/m 2 and 
around 1 nm in thickness indicates that the configura- 
tion described by Lubk et al. is very unlikely to occur in 
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FIG. 4. (Color online.) Same measures as in Fig. [3] but for 
the [111] (110) [111] 3/0 DW configuration (polarizations at an 
angle of 180°). 



reality. 

Among all the considered DW configurations, let us 
discuss in some detail two additional cases, namely, the 
[111](110) [111] 1/1 DW that is sketched in Fig.gc) and 
whose structural features are shown in Fig. [5j and the 
[111] (110) [111] 1/1 DW characterized by the results in 
Fig. |(3| Both of them are 71° FE-DWs, their respec- 
tive energies being 152 mJ/m 2 and 142 mJ/m 2 for the 
80-atom unit cell. Interestingly, while the former sat- 
isfies the m echanical compatibility condition mentioned 
abovep^the latter does not. Yet, in spite of this fact, 
the second DW has the lower energy. Indeed, generally 
speaking, our results show that the mechanical mismatch 
does not involve a large energy penalty. The relative 
small importance of this mechanical mismatch should 
probably be attributed to the small shear strains asso- 
ciated with the BFO's i?3c phase, which result in a small 
deformation of the DW plane of little energetic conse- 
quence. Hence, one should probably not consider the 
mechanical-matching criterion in studies of rhombohc- 
dral BFO, and should be cautious about its assumption 
for other materials or phases displaying small shears. 

From the three cases studied in detail so far, we see 
that low DW energies are associated to configurations in 
which the 0-0 bonds keep lengths similar to the ones 
displayed in bulk BFO. To further investigate this trend, 



9 



•< 0.25 
1 

tf) 

m "0-25 b 



•S. 0.25 -(b) 
I 

CO 

.? -0.25 - 



^ 10 

I 
d°-10 



10 

i °~ 



-< 



2.20 



2.00 ■ 

O 

<b 1.80 - 



-< 



0.20 



■g 0.10 — 
O 

O 0.00 1 



(a)" 



(d)_ 



(e) 



(f) 



8* 



* ft 



$ * » » 

A 



•I 



A 



If 



o_o_ 



0--O-0-0 



■l-J-ji-i-V-b- 



_x)--o--a-Q- 



o „ 



A 



A 



q_.q_.q__0. 

I I 



I 1 
* • * $ 

A-A-A__ A . 

I i 



I 1 

AA-A-A 
I i 



I 1 

I--H--H- 
a-H--H- 



I 1 

Q-Q-O-O-O --o- 

I I I I I 



X 


[100] 


+ 


[010] 


O 


[001] 



X 


[100] 


+ 


[010] 


o 


[001] 



o 


[100] 


* 


[010] 


A 


[001] 


■ i 







[100] 


# 


[010] 


A 


[001] 



-25 -20 -15 -10 -5 5 10 15 20 25 
Distance to DW center (A) 

FIG. 5. (Color online.) Same measures as in Fig. [3j but for 
the [111] (110) [111] 1/1 DW configuration (polarizations at an 
angle of around 71°). 
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FIG. 6. (Color online.) Same measures as in Fig. [3] but for 
the [ill] (110) [ill] 1/1 DW configuration (polarizations at an 
angle of around 71°; matching condition not fulfilled). 



we quantified the 0-0 distortion by measuring the dis- 
persion ri™ s of the N points in graphs such as that in 
Fig. pftf), relative to the values in bulk BFO, ^ ms , 



D l 



\ 



1 N 



4 ms ) 



(3) 



We computed this measure for each of the DW config- 
urations of Table [TXJ and plotted the DW energy as a 
function of it in Fig. [7J We see a very clear correlation 
between low DW energies and low variations in the 0-0 
bond lengths. 

We can further understand this trend by looking at the 
rotations that are present in the center of the domains, 
independently of the exact atomic rearrangements that 
take place at the DW. Let us consider for a moment a 
cubic perovskite crystal where we remove the A and B 
cations, and we keep the network of regular Oq octahe- 
dra. We now impose the formation of an AFD-DW, so 
that at each side of it the octahedra rotate a fixed angle, 
in a^a^a" fashion, about one of the (111) axes. The 
oxygen atoms at the DW plane are shared by two octa- 
hedra, one from each of the domains separated by the 
wall. Then, when the two different a~a~a~-like patterns 
freeze in, the rotation of each octahedra will try to send 



the oxygen atoms moving in one of four possible direc- 
tions that form 90° among them. Depending on which 
directions are chosen by the two octahedra sharing a sin- 
gle oxygen, three different cases can arise, as depicted in 
Fig. [7] the two directions may coincide (case A), they 
may form 90° (case B), or they may form 180° (case C). 
In order to accommodate the two conflicting movements 
required from a single oxygen ion, the octahedra at the 
DW will distort from their regular shape in the 90° case, 
and more so in the 180° case. 

In a material where AFD and FE distortions are 
present, the 06 octahedra are more distorted than in 
the case we have just described. However, if the rotation 
angle is relatively large, as in BFO, most of the displace- 
ment of the O atoms out of high-simmetry positions is 
due to the rotation about (111). We can therefore ap- 
ply the analysis of the previous paragraph to our BFO 
DWs, and label them following the same notation, since 
the kind of mismatch is the same for every pair of octa- 
hedra on opposite sides of the wall. We have done so in 
Table pi} As expected, we see that type A DWs are the 
lowest in energy, followed by type B DWs; finally, we find 
that type C DWs tend to be the most unfavorable ones. 
This is also illustrated in Fig. [7] 
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ductivity at room temperature in BFO DWs. In their ar- 
ticle, they also included results about the electronic prop- 
erties of the DWs obtained using first-principles calcula- 
tions; they mentioned that at the center of the domains 
the local density of states (LDoS) resembles the one of 
the bulk, while as the DW is approached the bandgap is 
reduced by an amount of up to 0.2 eV depending on the 
DW orientation. They correlated this reduction with the 
presence of the conductivity at the walls. 

Figure [8] shows the LDoS that we obtained from the 
atoms that are closest to the DW in each of the configu- 
rations of Tables IIIII and IIVI As we can see, if there are 



bandgap closings in these structures they are too small 
to be captured by this first-principles methodology. In 
these conditions, our interpretation is that conductivity 
at the DWs is not caused by any significative band-gap 
closing resulting from structural changes at the atomic 
level near a DW of pure BFO. 



FIG. 7. (Color online.) Left panel: DW energy as a function 
of a measure of distortions in O-O edges of 06 octahedra (see 
text); the symbols refer to different types of octahedra mis- 
matching at the DW. Right panel: Illustration of octahedra 
mismatches, as described in the text; each arrow indicates the 
axis around which each octahedron rotates. 



B. Role of exchange-correlation functionals 

In a previous article^ some of us pointed out how 
the use of different approximations to the exchange- 
correlation functional of DFT can affect the results of 
the total energy of BFO configurations. In particular, 
we found that the energy differences between phases (i.e., 
their relative stability) depends significantly on the em- 
ployed energy functional. Therefore, since we may in 
principle expect a similarly important dependence of the 
calculated DW energies, we repeated our calculations for 
three representative DW configurations usin g two other 
functionals: PBE 27 and PBEsolP Table Ev] shows small 
variations in the numbers obtained for the DW energies, 
while the relative ordering of the configurations is al- 
ways the same. The variations are smaller than those in 
Ref. [15] because our DW configurations are all relatively 
similar, while in the cited work phases with very dif- 
ferent structures (including, e.g., supertetragonal phases 
with O5 pyramids instead of Oq octahedra) were com- 
pared. We take these results as further confirmation that 
DFT predicts that the most stable DW configurations are 
those that keep the 0§ octahedra bond distortions to a 
minimum, and that their DW energies are in the range 
of several tens of mJ/m 2 . 



C. Electronic structure 

Seidel and collaborators^ were the first to report the 
experimental observation of an enhanced electronic con- 



D. Connection with experiment 

Seidel et a/P reported an enhanced conductivity in FE- 
DWs of types 2 (109°) and 3 (180°), while they found 
that FE-DWs of type 1 (71°) did not conduct in their 
experiments. On the other hand, the experiments of 
Farokhipoor and Noheda^ revealed an enhanced conduc- 
tivity at type 1 FE-DWs as well. The first set of authors 
argued that the observed conductivity was consistent 
with changes in the BFO structure at the domain wall; 
their complementary first-principles calculations showed 
a small reduction in the bandgap of pure BFO that they 
correlated to some degree with the conducting behavior. 
On the other hand, Farokhipoor and Noheda concluded 
that conduction at type 1 FE-DWs was related to the 
presence of oxygen vacancies. In another paper, Seidel 
and coworkers^ also demonstrated that the wall con- 
ductivity in La-doped BFO could be controlled through 
chemical doping with oxygen vacancies. 

Our results support a mechanism for the conductiv- 
ity that goes beyond the (small) effect that DW-specific 
structural distortions have in the local electronic proper- 
ties. As mentioned above, we find negligible differences 
between the electronic structure of the DWs and that of 
the bulk material. Further, we also investigated the ex- 
change couplings between iron spins at the DW, in order 



TABLE IV. Energies (in mJ/m 2 ) of the three representative 
DW configurations of Table |IH[ when studied with three dif- 
ferent approximations for the exchange-correlation functional 
(simulation cells with 80 atoms). 



DW Type LDA+f/ PBE+U PBEsol+t/ 



109° 


62 


78 


59 


180° 


73 


86 


91 


71° 


152 


136 


146 
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FIG. 8. (Color online.) LDoS of bulk BFO and of the 



three representative DW configurations of Table III (simu- 
lation cells with 80 atoms). The density of states has been 
projected on the closest 4 Fe, 4 Bi, and 12 O ions to the DW 
(thick lines), and on the d orbitals of those ions (shaded ar- 
eas). The bandgap is highlighted between two discontinuous 
vertical lines. 



to check whether a possible reduction in the magnitude of 
the magnetic interactions could have an effect on the spin 
structure (which might, e.g., be associated with a reduced 
Neel temperature at the DW) and thus on the electronic 
structure and conductivity; however, the observed vari- 
ations never exceeded a 10 % of the bulk exchange cou- 
plings, suggesting that the effect in the magnetic (and 
electronic) structure will be negligible. Hence, our cal- 
culations strongly indicate that the observed enhanced 
conductivity is to be linked to extrinsic factors such as 
oxigen vacancies or other sources of off-stoichiometry at 
the DWs. In this sense, we are currently performing cal- 
culations to compute the energy of formation of various 
defects at several types of DWs, wanting to determine (i) 
their effect in the electronic structure and (ii) their pref- 
erential occurrence in specific DWs. We believe that this 
is a promising route to explain the experimental obser- 
vations of a relatively large DW conductivity from first- 
principles theory. 

According to the results in Table [TTJ the creation of 
type 2 and type 3 DWs in bulk, defect-free BFO is favored 
energetically over the creation of type 1 DWs. On electri- 
cal switching at high field experiments, Seidel et aZPwere 



able to create the three types of domains on 100-nm-thick 
epitaxial BFO films grown on a (110) SrTi03 substrate. 
On the other hand, Farokhipoor and Noheda-' mention 
that their pulsed-laser deposition films grown on SrTiOs- 
(001) substrates showed mostly DWs of type 1. Such a 
preferential occurrence of 71° FE-DWs, which had also 
been observed in previous works is in obvious conflict 
with our preditions and points at the importance of fac- 
tors that were not included in our simulations. Indeed, 
the experimental results are for thin films grown in par- 
ticular orientations and conditions of misfit strain, and 
in the presence of a bottom substrate (and electrode) 
that plays a role as important as to induce a symme- 
try breaking in the possible orientations of the polar- 
ization (i.e., from the 8 (111) polarization variants that 
can in principle occur, the films of Ref. [5| only exhibit 4 
types of domains with the polarization always pointing 
towards the substrate). As regards the misfit strain, the 
recent first-principles investigation of Wei et al^ sug- 
gests that the DW hierarchy we obtained for the bulk 
case should be preserved for films grown on compressive 
SrTiO3-(001) substrates (with a misfit strain of about 
— 1.5%). Hence, the discrepancy between our predictions 
and the experimental observations seems to indicate that 
there are other factors beyond the epitaxial strain - e.g., 
the presence of defects, the effect of interfaces with sub- 
strate and/or electrodes, the specific electric boundary 
conditions - that play a key role in determining which 
DWs occur spontaneously in BFO films. 

Our first-principles results point to a picture of ex- 
tremely thin ferroelectric DWs around 1 nm thick, in 
agreement with previous calculations for other perovskite 
oxides. Our ideal DWs lie on infinite parallel planes, and 
they do not interact with each other when they are more 
than a few nanometers apart. (As evidenced by the rapid 
convergence of the computed DW energies as we increase 
the separation between DWs in our simulation box.) In- 
stead, real ferroelectrics show DWs in equilibrium that 
exhibit a characteristic roughness, S3 indicating that there 
are defects in the lattice that are pinning the wallsP In 
this respect, our simulations can be taken as the initial 
step in building a model of ferroelectric domains that 
would take into account the role of size effects and de- 
fects, a challenge that is currently out of reach of direct 
first-principles calculations. 



IV. SUMMARY AND CONCLUSIONS 

We have presented a theoretical study of domain walls 
(DWs) in perovskite oxides where two primary order pa- 
rameters coexist. We have taken BiFeOs as a sample case 
in which two structural instabilities of similar strength - 
i.e., a ferroelectric (FE) one, and one that involves con- 
certed rotations of the 06 octahedra (anti-ferrodistortive 
or AFD) - are present in the phase that is stable at am- 
bient conditions. We have used density-functional theory 
to perform atomic scale simulations of the hybrid domain 
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boundaries involving both FE and AFD discontinuities, 
and have thus revealed several novel aspects that need to 
be taken into account in situations like this one. First, 
and in contrast with what is implicitly assumed in most 
BFO studies, we have found that the DW energetics is 
essentially determined by the type of discontinuity in the 
AFD patterns, and not so much by the change in the 
FE polarization accross the wall. Thus, the AFD dis- 
tortions, which are usually treated as the secondary (i.e., 
less important) order parameter in BFO, turn out to play 
the leading role in determining which DWs are most sta- 
ble. Second, we have found that the lowest-energy DWs 
present atomic structures that can be directly identified 
with low-lying meta-stable BFO phases. Indeed, our re- 
sults suggest that multi-domain configurations in BFO 
are similar to polytypic structures with (meta-)stable 
polytypes occurring at the (very thin) domain bound- 
aries. This finding has an appealing and natural connec- 
tion with the novel nano-twinned structures^ that have 
been recently proposed to explain several controversial 
regions in the phase diagram of BFO and BFO-based 
solid solutions. Finally, we have discussed the implica- 
tions of our results as regards current experimental work 



on BFO's DWs, emphasizing the need of considering var- 
ious extrinsic factors to explain some of the most inter- 
esting experimental observations. 
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